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Abstract 

In recent contributions, algebraic multigrid methods have been designed and studied 
from the viewpoint of the spectral complementarity. In this note we focus our efforts 
on specific applications and, more precisely, on large linear systems arising from the 
approximation of weighted Laplacian with various boundary conditions. We adapt 
the multigrid idea to this specific setting and we present and critically discuss a wide 
numerical experimentation showing the potentiality of the considered approach. 



1 Introduction 



In the present note we test a specific application of a previously proposed 
algebraic multigrid procedure [?]. In that manuscript, we posed and partially 
answered the following question: having at our disposal an optimal multigrid 
procedure for A n x = b, {A n } being a given sequence of Hermitian positive 
definite matrices of increasing dimension, which are the minimal changes (if 
any) to the procedure for maintaining the optimality for B n y = c, {B n } new 
sequence of matrices with B n = A n + R n l 

Of course if there is no relation between {A n } and {B n } nothing can be said. 
However, under the mild assumption that there exists a value i? independent 
of n such that A n < "QB n and B n < MI n with M again independent of n, 
it has been clearly shown that the smoothers can be simply adapted and the 
prolongation and restriction operators can be substantially kept unchanged. 
The aim of this paper is to show the effectiveness of this approach in a specific 
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setting. More precisely, we consider linear systems A n (a)u = b arising from 
Finite Difference (FD) approximations of 

-V(o(z)Vu(s)) = f{x), x E tt = (0, l) d , d>l, 

where a(x) > clq > 0, f(x) are given bounded functions and with Dirichlet 
boundary conditions (BCs). Some remarks about the case of periodic or re- 
flective BCs are also considered (for a discussion on this topic see [?,?]). 
We recall that in the case a(x) = 1, the matrix A n (l) is structured, positive 
definite, ill-conditioned, and an optimal algebraic multigrid method is already 
available (see [?,?,?,?,?,?,?,?,?,?,?]) according to different BCs. 
Hereafter, owing to the spectral equivalence between the matrix sequences 
{A n (a}} and {A n (l)}, the key idea is that the multigrid procedure just de- 
vised for {A„(l)} can be successfully applied to {A n (a)} too. 
More in general in [?], we treated the case of structured-plus-banded uniformly 
bounded Hermitian positive definite linear systems, where the banded part R n 
which is added to the structured coefficient matrix A n is not necessarily defi- 
nite and not necessarily structured. In our setting A n = A n (l) is the structured 
part (it is Toeplitz, circulant etc, according to BCs) and R n = A n (a — 1) is 
the non-structured, non necessarily definite contribution. 
However, while a theoretical analysis of the Two-Grid Method (TGM) for 
structured+banded uniformly bounded Hermitian positive definite linear sys- 
tems has been given in [?], in terms of the algebraic multigrid theory by Ruge 
and Stiiben [?], the corresponding analysis for the multigrid method (MGM) 
is not complete and deserves further attention. Here, for MGM algorithm, we 
mean the simplest (and less expensive) version of the large family of multigrid 
methods, i.e., the V-cycle procedure: for a brief description of the TGM and 
of the V-cycle algorithms we refer to Section [2], while an extensive treatment 
can be found in [?], and especially in [?]. 

Indeed, the numerics in this note suggest that the MGM is optimal in the 
sense that (see [?]) the cost of solving the linear system (inverse problem) 
is proportional, by a pure constant not depending on n, to the cost of the 
matrix- vector product (direct problem): in our case more details can be given 
and in fact: 

a. the observed number of iterations is bounded by a constant independent 
of the size of the algebraic problem; 

b. the cost per iteration (in terms of arithmetic operations) is just linear 
as the size of the algebraic problem. 

Furthermore, given the spectral equivalence between {A n (a)}, a(x) > a > 0, 
and {t4„(1)}, a simpler numerical strategy could be used: use A n {l) as pre- 
conditioner for A n (a) in a PCG method and solve the linear systems with 
coefficient matrix A n (l) by MGM. Of course, this approach is simpler to im- 
plement, but since several linear systems have to be solved by MGM, the flop 
count can be more favorable in applying the MGM directly instead of using 
it as solver for the preconditioner. 
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The paper is organized as follows. In Section [2] we report the standard TGM 
and MGM algorithms, together with the reference theoretical results on the 
TGM optimal rate of convergence, under some general and weak assump- 
tions. In Section [3] the proposed approach is applied to the discrete weighted 
Laplacian and several numerical experiments are considered, by varying the 
diffusion function a(x) with respect to its analytical features. Finally, Section 
H] deals with further considerations concerning future work and perspectives. 



2 Two-grid and Multigrid Method 



We carefully report the TGM and MGM algorithms and we describe the the- 
oretical ground on which we base our proposal. We start with the simpler 
TGM and then we describe the MGM and its interpretation as stationary or 
multi-iterative method, see [?]. 



2.1 Algorithm definition 



Let n be a positive ci-index, d > 1, and let N(-) be an increasing function 
with respect to no- In devising a TGM and a MGM for the linear system 
A no x no = b no , where A no G C N{no)xN{no) and x no ,b no E C N{n °\ the ingredients 
below must be considered. 

Let ni < no (componentwise) and let p 7 ^ £ £%o)xJV(ni) a gj ven full-rank 
matrix. In order to simplify the notation, in the following we will refer to any 
multi-index n s by means of its subscript s, so that, e.g. A s := A ng , b s := b ns , 
P s s +1 '=1%? 1 , etc. 

With these notations, a class of stationary iterative methods of the form 
x (j+i) = V s x^p +b s is also considered in such a way that Smooth(x^\ b s , V s , u s ) 
denotes the application of this rule v s times, with v s positive integer number, 
at the dimension corresponding to the index s. 

Thus, the solution of the linear system A no x no = b no is obtained by applying 
repeatedly the TGM iteration, where the j th iteration 

x =TGM(x ( J> ,b ,A Q ,V 

,pre> ^Ojpost; ^O.post) 
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is defined by the following algorithm [?]: 

y := TGM(x , b , A , V ,pre; ^O.pre; Vb,post; ^O.post) 



x := Smooth(x , b , V QyPie , z/ , P rc) 



Pre- smoothing iterations 



r := b - A x 
n := {pl) H ro 

Solve A x y x = n, with A x := (pl) H A Po 



Vo 



xo + PoVi 



Exact Coarse Grid Correction 



yo := Smooth(y , b , V 0:POSt , ^o.post) 



Post-smoothing iterations 



The first and last steps concern the application of u 0tPrc steps of the pre- 
smoothing (or intermediate) iteration and of ^o )POS t steps of the post- smoothing 
iteration, respectively. Moreover, the intermediate steps define the so called 
coarse grid correction, that depends on the projection operator {po) H ■ In such 
a way, the TGM iteration represents a classical stationary iterative method 
whose iteration matrix is given by 



TGM 



V'-^Ojpost f^r^f^ T/"^0,prc 
0,post <~^<^0 ^0,pre ? 



(2-1) 



where CGC = h - pi (pl) H 'AqpI 



(pl) H A denotes the coarse grid correc- 
tion iteration matrix. 

The names intermediate and smoothing iteration used above refer to the multi- 
iterative terminology [?]: we say that a method is multi- iterative if it is com- 
posed by at least two distinct iterations. The idea is that these basic compo- 
nents should have complementary spectral behaviors so that the whole pro- 
cedure is quickly convergent (for details see [?] and Sections 7.2 and 7.3 in 
[?]). Notice that in the setting of Hermitian positive definite and uniformly 
bounded sequences, the subspace where A is ill-conditioned corresponds to 
the subspace in which A has small eigenvalues. 

Starting from the TGM, the MGM can be introduced as follows: instead of 
solving directly the linear system with coefficient matrix A\, the projection 
strategy is recursively applied, so obtaining a multigrid method. 
Let us use the Galerkin formulation and let n > n\ > . . . > rii > 0, with I 
being the maximal number of recursive calls and with N(n s ) being the corre- 
sponding matrix sizes. 

The corresponding MGM generates the j th iteration 



Xq + -* — M.GM(0, Xq \ b , A , Vo )Pre , ^0,pre, K),post) ^0,posty 
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according to the following algorithm: 

y s . AAGM^ (s, X s , b s , A SJ V^p re , i^pre; ^fjpost; ^s,post) 



if s = I then 



Solve(A s y s = b s 



Exact solution 



else 



X s := Smooth (x s , b s , K.pre, ^s.pre) 



Pre- smoothing iterations 



r s := b s — A s x s Coarse Grid Correction 

r s+i ■= (p s s +1 ) H r s 
y s +i-= MGM(s + 1,0 

s+li ^s+lj^s+l.prej ^s+ljpre; ^s+l.post; ^s+ljpost) 

y s := x s +p s s +1 y s+1 



y s \= Smooth(y s ,b s , 



V, 



s,post; ^s,post 



Post- smoothing iterations 



-1\H 



A s p s s +1 is more profitably computed in the 



where the matrix A s+ \ := (p 
so called pre- computing phase. 

Since the MGM is again a linear fixed-point method, the j th iteration Xq^ ^ 
can be expressed as MGMqXq^ +(/o — MGMq)Aq 1 bo, where the iteration 
matrix MGMq is recursively defined according to the following rule (see [?]): 



MGMi 
MGM, 



O, 

T 7"Ks .post 

^s,post 



Is-p s s + Vs + i-MGM s+1 ]A;IM +1 ) h A 

s = 0, 



T/"Ks ,prc 

*s,pre 5 

..J-h 



(2.2) 



and with MGM S and MGM s+ i denoting the iteration matrices of the multi- 
grid procedures at two subsequent levels. 

At the last recursion level Z, the linear system is solved by a direct method 
and hence it can be interpreted as an iterative method converging in a single 
step: this motivates the chosen initial condition MGM\ = O. 
By comparing the TGM and MGM, we observe that the coarse grid correc- 
tion operator CGC S is replaced by an approximation, since the matrix A~^_± 
is approximated by (I s +i -MGM S+1 ) A~^_ x as implicitly described in (12. 2p for 
s = 0, ...,/ — 1. In this way step 4., at the highest level s = 0, represents an 
approximation of the exact solution of step 4. displayed in the TGM algorithm 
(for the matrix analog compare (12. 2p and (12.11) ). Finally, for / = 1 the MGM 
reduces to the TGM if Solve(Aiy\ = b\) is y\ = A[ b%. 
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2.2 Some theoretical results on TGM convergence and optimality 



In this paper we refer to the multigrid solution of special linear systems of the 
form 

B n x = b, B n E C N{n)xN{n \ x, b G C N{n) (2.3) 

with {B n } Hermitian positive definite uniformly bounded matrix sequence, n 
being a positive d-index, d > 1 and N(-) an increasing function with respect 
to it. More precisely, we assume that there exists {A n } Hermitian positive 
definite matrix sequence such that some order relation is linking {A n } and 
{£>„}, for n large enough and we suppose that an optimal algebraic multigrid 
method is available for the solution of the systems 

A n x = 6, A n G C N{n)xN{n \ x, b G C N(n \ (2.4) 

The underlying idea is to apply for the systems (12.31) the some algebraic TGM 
and MGM considered for the systems (12.41) . i.e., when considering the very 
same projectors. In fact, the quoted choice will give rise to a relevant simpli- 
fication, since it is well-known that a very crucial role in MGM is played by 
the choice of projector operator. 

In the algebraic multigrid theory some relevant convergence results are due to 
Ruge and Stiiben [?], to which we referred in order to prove our convergence 
results. 

Hereafter, by || • || 2 we denote the Euclidean norm on C m and the associated 
induced matrix norm over C mxm . If X is Hermitian positive definite, then its 
square root obtained via the Schur decomposition is well defined and positive 
definite. As a consequence we can set || • \\x = WX 1 ^ 2 ■ W2 the Euclidean norm 
weighted by X on C m , and the associated induced matrix norm. In addition, 
the notation X < Y, with X and Y Hermitian matrices, means that Y — X 
is nonnegative definite. In addition the sequence {A n }, with X n Hermitian 
positive definite matrices, is a uniformly bounded matrix sequence if there 
exists M > independent of n such that ||X n ||2 < M, for n large enough. 

Theorem 2.1 [?] Let Aq be a Hermitian positive definite matrix of size N(no), 
let pi G C Ar ( n °) x7V ( ni ) ; no > ni, be a given full-rank matrix and let Vo. pos t be 
the post- smoothing iteration matrix. Suppose that there exists a post > ; inde- 
pendent of no, such that for all x G C^^ ^ 

ll^O.post^H^o < \\ x \\aq ~ "post IMIaoD^Ao' ( 2 ^) 

where Dq is the diagonal matrix formed by the diagonal entries of Aq. 
Assume, also, that there exists (3 > 0, independent of n^, such that for all 
x G C N ^ 

y ™&J x -poy\\°°^P II*- ( 2 - 6 ) 



Then, (3 > a post and \\TGM \\ Ao < Jl - a post //3 < 1 
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Notice that all the constants a pos t and (3 are required to be independent of 
the actual dimension in order to ensure a TGM convergence rate independent 
of the size of the algebraic problem. 

It is worth stressing that Theorem 12.11 still holds if the diagonal matrix D 
is replaced by any Hermitian positive matrix X Q (see e.g. [?]). Thus, X — I 
could be a proper choice for its simplicity. 

Thus, by referring to the problem in 12.31 we can claim the following results. 

Proposition 2.2 [?] Let {A n } be a matrix sequence with A n Hermitian posi- 
tive definite matrices and let p^ E (£N{no)xN(m) ^ e a gi ven full-rank matrix for 
any n Q > such that there exists (3a > independent of n so that for all 
x E C N ^ 

mm \\x - ply\\l < (3 a \\x\\a - ( 2 -7) 

Let {B n } be another matrix sequence, with B n Hermitian positive definite ma- 
trices, such that A n < "8B n , for n large enough, with d > absolute constant. 
Then, for all x E C^™ -* and n large enough, it also holds (3b = Pa$ and 

min , \\x - plvWl < Pb\\x\\1 - (2.8) 



Therefore, the convergence result in Theorem 12.11 holds true also for the matrix 
sequence {B n }, if the validity of condition (12.51) it is also guaranteed. It is worth 
stressing that in the case of Richardson smoothers such topic is not related to 
any partial ordering relation connecting the Hermitian matrix sequences {A n } 
and {Bn}, i.e. inequalities (12.51) . and the corresponding for the pre-smoother 
case, with {-£>„} instead of {A„}, have to be proved independently. 

Proposition 2.3 [?] Let {B n } be an uniformly bounded matrix sequence, 
with B n Hermitian positive definite matrices. For any n > 0, let V riyWe = 
I n — u pre B n , V^ )POSt = I n — u post B n be the pre-smoothing and post-smoothing 
iteration matrices, respectively considered in the TGM algorithm. Then, there 
exist ensure, «s,post > independent of such that for all x E C^ 1 -™ ^ 

\\Vo (2.9) 
\\Vo iPOS tx\\ Bo < \\x\\ 2 Bo - a B>post \\x\\ 2 B 2. (2.10) 

See Proposition 3 in [?] for the analogous claim in the case of z/ pre , v post > 0. 

In this way, according to the Ruge and Stiiben algebraic theory, we have proved 
the TGM optimality, that is its convergence rate independent of the size N(n) 
of the involved algebraic problem. 

Theorem 2.4 [?] Let {B n } be an uniformly bounded matrix sequence, with B n 
Hermitian positive definite matrices. Under the same assumptions of Proposi- 
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Hons 12.21 and 12.31 the TGM with only one step of post-smoothing converges to 
the solution of B n x = b and its convergence rate is independent of N(n). 

Clearly, as just discussed in [?], the TGM iteration with both pre-smoothing 
and post-smoothing is never worse than the TGM iteration with only post- 
smoothing. Therefore Theorem 12.41 implies that the TGM with both post- 
smoothing and pre-smoothing has a convergence rate independent of the di- 
mension for systems with matrices B n under the same assumptions as in The- 
orem [231 

Furthermore, the same issues as before, but in connection with the MGM, 
deserve to be discussed. First of all, we expect that a more severe assumption 
between {A n } and {B n } has to be fulfilled in order to infer the MGM opti- 
mality for {B n } starting from the MGM optimality for {A n }. The reason is 
that the TGM is just a special instance of the MGM when setting / = 1. 
In the TGM setting we have assumed a one side ordering relation: here the 
most natural step is to consider a two side ordering relation, that is to as- 
sume that there exist positive constants i?i,i?2 independent of n such that 
$iB n < A n < $2-Bn> for every n large enough. The above relationships simply 
represent the spectral equivalence condition for sequences of Hermitian pos- 
itive definite matrices, which is plainly fulfilled in our setting whenever the 
weight function is positive, well separated from zero, and bounded. 
In the context of the preconditioned conjugate gradient method (see [?]), it 
is well known that if {P n } is a given sequence of optimal (i.e., spectrally 
equivalent) preconditioners for then {P n } is also a sequence of optimal 

preconditioners for {B n } (see e.g. [?]). The latter fact just follows from the 
observation that the spectral equivalence is an equivalence relation and hence 
is transitive. 

In summary, we have enough heuristic motivations in order to conjecture that 
the spectral equivalence is the correct, sufficient assumption and, in reality, 
the numerical experiments reported in Section [3] give a support to the latter 
statement. Refer to [?] for some further remark about this topic. 



3 Numerical Examples 

Hereafter, the aim relies in testing our TGM and MGM (standard V-cycle 
according to Section [2]) applied to standard FD approximations to 

- V(a(x)Vw(x)) = f(x), x E tt = (0, d > 1, (3.1) 

with assigned BCs and for several choices examples of the diffusion coefficient 
a(x) > ao > 0. 

The projectors are properly chosen according to the nature of structured part, 
that depends on the imposed BCs. For instance, in the case of Dirichlet BCs 
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we split the arising FD matrix A n (a) as 

A n (a) = a min T n (A n (l)) + R n (a), R n {a) = A n (a) - a min T n (A n (l)) , 

where r n (A n (l)) denotes the FD matrix belonging to the r (or DST-I) algebra 
[?] obtained in the case of a(x) = 1 and a min equals the minimum of a(x) on 
Cl in order to guarantee the positivity of R n (a). 

On the other hand, we will use, in general as first choice, the Richardson 
smoothing/intermediate iteration step twice in each iteration, before and af- 
ter the coarse grid correction, with different values of the parameter uj. In some 
cases better results are obtained by considering the Gauss-Seidel method for 
the pre-smoothing iteration. 

According to the algorithm in Section [21 when considering the TGM, the exact 
solution of the system is obtained by using a direct solver in the immediately 
subsequent coarse grid dimension, while, when considering the MGM, the ex- 
act solution of the system is computed by the same direct solver, when the 
coarse grid dimension equals 15 d (where d = 1 for the one-level case and d = 2 
for the two- level case). 

In all tables we report the numbers of iterations required for the TGM or 
MGM convergence, assumed to be reached when the Euclidean norm of the 
relative residual becomes less than 10~ 7 . We point out that the CPU times 
are consistent with the iteration counts. 

Finally, we stress that at every level (except for the coarsest) the structured 
matrix parts are never formed since we need only to store the nonzero Fourier 
coefficients of the generating function at every level for matrix-vector multi- 
plications. Thus, besides the 0(N(n)) operations complexity of the proposed 
MGM both with respect to the structured part and clearly with respect to the 
non-structured one, the memory requirements of the structured part are also 
very low since there are only 0(1) nonzero Fourier coefficients of the generat- 
ing function at every level. On the other hand, the projections of the initial 
matrix correction R n (a) are stored at each level according to standard sparse 
matrix techniques during the pre-computing phase. 

3.1 DirichletBCs 

We begin by considering the FD approximation of (13.11) with Dirichlet BCs 
in the one-level setting. As already outlined, in this case the arising matrix 
sequence {A n (a)} can be split as 

A n (a) = a min r n (A n (l)) + R n (a), R n (a) = A n (a) - a min r n (A n (l)), 

where T n (A n {l)) and a m i n are defined as before. More precisely, {r n {A n {l))} is 
the T/Toeplitz matrix sequence generated by the function f{t) = 2 — 2 cos(t), 
t G (0, 2tt] and a min equals the minimum of a(x) on Cl. 
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Let us consider Aq(o) G IR n ° xn °, with 1-index no > (according to the notation 
introduced in Section[2j we refer to any multi-index n s by means of its subscript 
s). Following [?,?], we denote by Tq G fl^o*"^ n ^ _ _|_ ^ ^he operator such 
that 

1 for i — 2j, ?' = 1,..., m, 
(tf)ij={ (3.2) 
I otherwise, 

and we define a projector (pl) H , pi G ]R n ° xni as 

= 4^oT \ P = tridiag [1, 2, 1] = r (/), f(t) =2 + 2 cos(t). (3.3) 

On the other hand, for the smoothing/intermediate Richardson iterations, the 
parameters u are chosen as 

k'pre 2/(||/||oo+||Pn(a)||oo) 
^post=l/(||/||oo+||i?n(a)||oo), 

and we set u prc = z/ post = 1. 

The first set of numerical tests refer to the following settings: a(x) = 1, 
a(x) = e x , a(x) = e x + 1, (denoted in short as al, a2, a3 respectively). 
In Table [1] we report the numbers of iterations required for the TGM con- 
vergence, both in the case of the Richardson pair, and of the Richardson + 
Gauss-Seidel pair. All these results confirm the optimality of the proposed 
TGM in the sense that the number of iterations is uniformly bounded by a 
constant not depending on the size N(n) indicated in the first column. 
In Table [2] we report the some results, but with respect to the V-cycle applica- 
tion. The numerics seems allow to claim the optimality convergence property 
can be extended to the MGM. 

It is worth stressing that the difference in considering Richardson or Gauss- 
Seidel in the pre-smoothing iterations is quite negligible in the MGM case. 
In Table |3] we report a deeper analysis of the TGM superlinear behavior in 
the a2 setting. More precisely, we consider the test functions a(x) = e x + 10 k 
with k ranging from to 6. The convergence behavior is unaltered in the case 
of the Richardson + Gauss-Seidel pair, while for increasing k we observe that 
the number of required iterations by considering the Richardson+Richardson 
pair progressively approaches the reference al case. In fact as k — >• oo the 
function a(x) after a proper scaling converges to the constant 1. 

The projector definition plainly extends to the two-level setting by using 
tensor arguments: {j>\) H is constructed in such a way that 

Po = PoUl (3.4) 
P = tridiagm [1,2,1]® tridiagm [1,2,1], (3.5) 

"o 

U^T^^T^) (3.6) 
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Table 1 

Number of iterations required by TGM - one-level case with Dirichlet BCs 



Richardson+ Richardson 




Richardson+Gauss- 


Seidel 


N(n) 


al 


a2 


a3 




N(n) 


al 


a2 


a3 


31 


2 


8 


5 




31 


8 


8 


8 


63 


2 


6 


4 




63 


8 


8 


8 


127 


2 


5 


4 




127 


8 


8 


8 


255 


2 


4 


4 




255 


8 


8 


8 


511 


2 


4 


3 




511 


8 


8 


8 



Table 2 

Number of iterations required by MGM - one-level case with Dirichlet BCs 



Richardson-I- Richardson 




Richardson+Gauss-Seidel 


N(n) 


al 


a2 


a3 


N(n) 


al 


a2 


a3 


15 


1 


1 


1 




15 


1 


1 


1 


31 


2 


8 


5 




31 


8 


8 


8 


63 


7 


7 


7 




63 


9 


9 


9 


127 


8 


8 


8 




127 


9 


9 


9 


255 


8 


8 


8 




255 


9 


9 


9 


511 


8 


8 


8 




511 


9 


9 


9 



Table 3 

Number of iterations required by TGM - one-level case with Dirichlet BCs 



Richardson-I- Richardson 




Richardson+Gauss 


-Seidel 






a(x) 


= e x + 10 fe 




a(x) 


= e x + 10 fe 


k 


k 


N(n) 


al 





1 


2 


3 


4 


5 


N(n) 


al 





1 


2 


3 


4 


5 


31 


2 


5 


4 


3 


3 


3 


2 




31 


8 


8 


8 


8 


8 


8 


8 


63 


2 


4 


4 


3 


3 


3 


2 




63 


8 


8 


8 


8 


8 


8 


8 


127 


2 


4 


4 


3 


3 


3 


2 




127 


8 


8 


8 


8 


8 


8 


8 


255 


2 


4 


3 


3 


3 


3 


2 




255 


8 


8 


8 


8 


8 


8 


8 


511 


2 


3 


3 


3 


3 


2 


2 




511 


8 


8 


8 


8 


8 


8 


8 



with 4 } = InV + 1 and where T \n^) G is the one-level matrix 

given in ( 13. 21) . 

The quoted choice represents the most trivial extension of the one-level pro- 
jector to the two-level setting and is also the less expensive from a compu- 
tational point of view: in fact, pi = t ((2 + 2cos(ii)(2 + 2 cos^)))^ equals 
[r n(1) (p(2 + 2cos(t 1 )))T 1 (4 1) )] ® [r nC2) (p(2 + 2 cos(t 2 )))T^n^)}. 

"o 

Tables H] and [5] report the number of iterations with the same notation as 
before and where we are considering the following function tests: a(x) = 1, 
a(x) = e xl+X2 , a(x) = e xl+X2 + 2, (denoted in short as al, a2, a3, respectively). 
Though the convergence behavior in the case of the Richardson+Richadson 
pair is quite slow, we can observe that the number of MGM iterations re- 
quired to achieve the convergence is essentially the same as in the TGM. This 
phenomenon is probably due to some inefficiency in considering the approx- 
imation 1 1 i?n (a) || oo i n t ne tuning of the parameter u pre and uj post . In fact, it 
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Table 4 

Number of iterations required by TGM - two-level case with Dirichlet BCs 



Richardson+ Richardson 




Richardson+Gauss-Seidel 


N(n) 


al 


a2 


a3 


N(n) 


al 


a2 


a3 


31 2 


16 


73 


38 




31 2 


13 


14 


14 


63 2 


16 


82 


41 




63 2 


13 


15 


14 


127 2 


16 


86 


43 




127 2 


13 


15 


14 


255 2 


16 


89 


44 




255 2 


13 


15 


14 



Table 5 

Number of iterations required by MGM - two-level case with Dirichlet BCs 



Richardson-I- Richardson 




Richardson+Gauss-Seidel 


N(n) 


al 


a2 


a3 


N(n) 


al 


a2 


a3 


15 2 


1 


1 


1 




15 2 


1 


1 


1 


31 2 


16 


73 


38 




31 2 


13 


14 


14 


63 2 


16 


83 


42 




63 2 


13 


15 


15 


127 2 


16 


88 


43 




127 2 


13 


15 


15 


255 2 


16 


90 


44 




255 2 


13 


15 


15 



is enough to substitute, for instance, the pre-smoother with the Gauss-Seidel 
method in order to preserve the optimality both in the TGM and the MGM 
case. 

Finally, in Table [6j we report the number of iterations required by MGM, 
in the case of some other test functions. More precisely, we are considering 
the C 1 function a(x,y) = e x+ \ y - l,2 ^'\ the C° function a(x,y) = e x+l ^ 1/21 , 
and the piecewise constant function a(x,y) — 1 if x, y < 1/2, S otherwise, 
with 5 = 10, 100, 1000 (denoted in short as a4, a5, a6, a7, and a8, respec- 
tively). Taking into account the previous remarks, our smoothing choice is 
represented by the Richardson+Gauss-Seidel pair. Moreover, the CG choice 
is also investigated, both in connection to the Richardson or the Gauss-Seidel 
smoother. 

The MGM optimality is again observed, according to a proper choice of the 
smoother pair. 

In conclusion for keeping a proper optimal convergence, we can claim that 
Gauss-Seidel is necessary and the best pair is with conjugate gradient. The 
explanation of this behavior is again possible in terms of multi-iterative pro- 
cedures and spectral complementarity: in fact while Richardson is effective 
essentially only in the high frequencies space, both Gauss-Seidel and CG are 
able to reduce the error also in the middle frequencies and in addition they 
are robust with respect to the scaling produced by the weight function a. 

3.2 Periodic and Reflective BCs 

Hereafter, we briefly address the case of periodic or reflective BCs. In par- 
ticular we focus on the structured part of the splitting related to the FD 
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Table 6 

Number of iterations required by MGM - two-level case with Dirichlet BCs (f 
more than N(n) iterations required for convergence) 



Richardson+Gauss-Seidel 



N(n) 


a4 


a5 


a6 


a7 


a8 


15 2 


1 


1 


1 


1 


1 


31 2 


14 


14 


13 


13 


13 


63 2 


15 


15 


13 


13 


13 


127 2 


15 


15 


14 


14 


14 


255 2 


15 


15 


14 


14 


14 



Richardson+ C G 


N(n) 


a4 


a5 


a6 


a7 


a8 


15 2 


1 


1 


1 


1 


1 


31 2 


21 


24 


46 


1472 


t 


63 2 


26 


28 


59 


1990 


t 


127 2 


26 


30 


64 


1783 


t 


255 2 


27 


31 


60 


1973 


t 



Gauss-Seidel+CG 



N(n) 


a4 


a5 


a6 


a7 


a8 


15 2 


1 


1 


1 


1 


1 


31 2 


12 


12 


11 


10 


10 


63 2 


12 


12 


11 


10 


10 


127 2 


12 


12 


11 


10 


10 


255 2 


12 


12 


11 


10 


10 



discretization with respect to a(x) = 1, since our multigrid strategy is tuned 
just with respect to it. 

In the case of periodic BCs the obtained matrix sequence is the one-level cir- 
culant matrix sequence {S n (f)} generated by the function f(t) = 2 — 2 cos(t), 



t G (0, 2tt]. Following [?], we consider the operator Tq e 
such that 

1 for i = 2j - 1, j = 1, 



tmoxrii 



n 



2ni, 



(To 1 ; 



otherwise, 

and we define a projector (pl) H , pi G R™° xril , as pi = PqTq, P = So(p), pit) = 
2 + 2cos(t). Clearly, the arising matrices are singular, so that we consider, for 
instance, the classical Strang correction [?] 

SnM) = SnM) + i ' (a^ot) ^wr 

where e is the vector of all ones. 

By using tensor arguments, our approach plainly extend to the two-level set- 
ting. 

When dealing with reflective BCs, the obtained matrix sequence is the one- 
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level DCT III matrix sequence C n (f) n generated by the function f(t) = 
2 - 2cos(t), t E (0,2?r]. Following [?], we consider the operator e R n ° xn \ 
riQ = 2ni, such that 



and we define a projector (p^) H , p l e R n ° xni , as p\ = P T^, P = C (p), 
p(t) = 2 + 2cos(£). Clearly, due to the singularity, we consider, for instance, 



Again, the two-level setting is treated by using tensor arguments. 
The numerical tests performed in the case of periodic or reflective BCs have 
the same flavor as those previously reported in the case of Dirichlet BCs and 
hence we do not report them since the observed numerical behavior gives the 
same information as in the case of Dirichlet BCs. 



4 Concluding Remarks 

We have presented a wide numerical experimentation concerning a multigrid 
technique for the discrete weighted Laplacian with various BCs . In accordance 
with the theoretical study in [?], the choice of the smoothers can be done 
taking into account the spectral complementarity, typical of any multi-iterative 
procedure. In particular, we have noticed that when the weight function a adds 
further difficulties in the middle frequencies (e.g., when a is discontinuous), 
the use of pure smoothers like Richardson, reducing the error only the high 
frequencies, is not sufficient. Conversely, both CG and Gauss-Seidel work also 
reasonably well in the middle frequencies (what is called the intermediate space 
in a multi-iterative method) and in fact, in some cases, their use is mandatory 
if we want to keep the optimality of the method, i.e., a convergence within a 
given accuracy and within a number of iterations not depending on the size 
of the considered algebraic problem. 




for i e {2j - l,2j}, j = l,...,ni, 
otherwise, 




14 



